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ABSTRACT 

Present numerical models describing pollutant diffusion 
employ two main approaches to the problem. Gaussian bivar- 
iate distributions based on the statistical properties of 
atmospheric flow have the advantages of simplicity and com- 
putational economy. Numerical simulations of diffusion, 
based on the continuity equation and conservation relations, 
allow a more precise description of atmospheric turbulent 
flow. 

The thesis examines numerical models of each method. 

The governing equations are discussed, with basic assumptions 
which permit reasonable simplifications to be made. Results 
of both models are examined qualitatively to describe pat- 
terns of marine-layer pollutant dispersal. Calculations 
from both models reveal wind velocity to be the most criti- 
cal atmospheric characteristic controlling dispersion. 

Rates of spreading are also closely related to other measur- 
able meteorological properties such as thermal stability and 
mixing depth. 
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I. INTRODUCTION 



The use of computer simulation techniques has in recent 
years greatly enhanced the exploitation of existing theoret- 
ical knowledge in the analysis of physical phenomena. The 
development and application of numerical predictive models 
is now routine in many scientific and technical fields, in- 
cluding airshed management. Such models can be used to study 
the complicated relationship between air quality and emis- 
sion sources as a function of meteorological parameters, 
topography and time. 

The goal of research in this field is the accurate pre- 
diction of pollution levels and, ultimately, air quality 
management. However, before prevention of harmful pollution 
episodes by regulation of sources is feasible, an adequate 
real-time forecast system must be developed. In addition to 
"red-flagging" pollutant concentrations in excess of levels 
toxic to living organisms, this scheme should have the ca- 
pacity to designate: 

1. the optimum loci and times of day at which atmos- 
pheric injection of pollutants will be least damaging to 
air quality; 

2. the relative effectiveness of various control meas- 
ures on a particular source. 

The success of such a control model on either a local or 
regional basis depends primarily upon the reliability of 
diffusion estimates calculated by the method- ’ Given the 
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current state of diffusion modeling, the achievement of these 
goals is not likely for some time to come. Present models 
are formulated from the concentration equation governing the 
various pollutant species. In theory this equation can be 
solved using known source emission rates, meteorological 
parameters, horizontal and vertical boundary conditions, 
turbulent transport rates and chemical reactivity. Practi- 
cally, however, present data acquisition falls far short of 
the necessary detail for a valid solution to the equation. 

As a consequence, various restricted but more tractable 
models have been developed, based upon statistical theory 
and simplifying assumptions. The two most common approaches 
today are Gaussian diffusion and numerical simulation. 

The Gaussian plume (steady-state) and puff (dynamic) 
models describe the concentration distribution of inert con- 
stituents and are derived from the stochastic continuity 
equation under homogeneous, isotropic conditions. Down- 
stream from a point source the average concentration on a 
crosswind plane is assumed to fit a normal bivariate statis- 
tical distribution (Figure 1) . Observations indicate that 
this approximation is reasonable only for very short-term 
forecasts. However, the shortage of detailed meteorological 
data and the simplicity and computational economy of Gaus- 
sian methods have motivated their continued popularity. The 
models of Pasquill [1961], Turner [1964], Hilst [1967], and 
Johnson [1970] are representative examples. 
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FIGURE 1 GAUSSIAN PLUME MODEL 
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Most statistical formulae are limited in their applica- 
tion by failure to consider chemical alterations of pollu- 
tant species. Some atmospheric constituents react chemical- 
ly or photochemically with others, yielding entirely new com- 
pounds. Deposition and absorption at the earth's surface 
and washout by aerosols reduce concentrations. Further, the 
normal statistical distribution accurately applies only to 
continuous sources emitting into a steady laminar flow. 
Conditions of zero wind or weak winds can result in singu- 
larities. The complex nature of thermally- or topograph- 
ically-induced turbulence and wind shear can only be crudely 
approximated for incorporation in Gaussian models. 

Atmospheric pollutant dispersion is more adequately de- 
scribed by numerical simulation of concentration diffusion. 
Simplified forms of the continuity equation have been applied 
by Randerson [1970], Lamb [1971], and Shir et al [1973] to 
form the bases for multiple source urban diffusion models. 

As yet, however, the closure problem encountered in limiting 
the degrees of freedom of the turbulent system still looms 
as an insurmountable problem. 
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II. PURPOSE 



The primary objective of this thesis is to study the im- 
portance of meteorological conditions in the prediction of 
pollutant dispersal. Variations of wind velocity, atmos- 
pheric stability and mixing-layer thickness are considered. 
Meteorological effects on dispersion are examined for both 
a simple Gaussian formula and for a model based on a more 
complete form of the concentration diffusion equation. 

Most present models distinguish between area and line 
sources of pollution, such as clustered home heating units 
and freeways, from the emissions of major atmospheric out- 
falls such as power plant stacks. Area and line contribu- 
tions can be regarded as the integral of numerous small 
point sources over a particular region and as such are sep- 
arable from large-point contributions. 

The single most significant cause of air pollution in 
many cities is the generation of electrical power by combus- 
tion of fossil fuels. This thesis therefore considers con- 
taminants from a large point source only, neglecting area 
and line contributions. Furthermore, since nearly fifty 
percent of all power produced in the United States is gen- 
erated on the coastal plain, the study is applied to that 
zone (specifically the Monterey Bay region) . 
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III. DESCRIPTION OF THE COASTAL ZONE 



In the transition zone between marine and continental 
airsheds, both natural inhomogeneities in temperature, den- 
sity and moisture, and human activities contribute to the 
concentration of atmospheric contaminants. An airmass flow- 
ing across the coastline encounters a step change in surface 
roughness and frequently in temperature as well. The dynam- 
ical and thermal response of the atmosphere to these changes 
in the properties of the underlying terrain is the down- 
stream propagation of turbulence in the boundary layer. 

The rate of diffusion of airborne pollutants, particu- 
larly in the vertical, depends primarily upon the turbulent 
character of the atmosphere. Turbulence is controlled by 
the wind velocity, static stability, nature of terrain, de- 
gree of surface heating, lateral boundaries, existence of 
inversion layers and other factors of local significance. 

The result is a reduction in mixing time in the coastal zone 
marine layer, enhancing the distribution of contaminants. 

A. MONTEREY BAY CIRCULATION 

The geography and topographic relief of the Monterey Bay 
region are shown in Figure 2. Although this study is theo- 
retical in that no attempt is made to correlate findings 
with actual source inventory data or measured pollution 
levels, the thesis models are designed to represent as 
closely as possible the meteorological features of the locale. 
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FIGURE 2 



SOUTH 



MONTEREY 
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Darling [1971] has shown that for the South Monterey Bay 
the low-level atmospheric circulation is almost totally to- 
pography-dependent, with little contribution by raeso- and 
synoptic-scale circulatory patterns. Read [1972] found the 
surface wind field at Moss Landing to be essentially repre- 
sentative of the undisturbed offshore flow. Circulation in 
the north and south parts of the bay is dominated by diver- 
gent topographically-induced eddies. 

The characteristic daytime thermal structure of the sur- 
face boundary layer exhibits a dry inversion near 600 meters 
throughout much of the year. Cool marine air intrudes along 
the coastline in the early morning, wedging the faster- 
warming continental airmass upward. The inversion forms at 
the interface between these two bodies of air and is usually 
intensified by subsidence of warm dry air from the eastern 
cell of the Pacific Subtropical Ridge. The turbulent, well- 
mixed marine layer increases in depth early in the day, 
reaching a maximum thickness first at the coastline and later 
inland; depth of the layer decreases in the afternoon. 
Schroeder et al [1967] related the depth of the marine air- 
mass directly to the wind velocity. 

The predominant local wind and thermal patterns thus af- 
fect the dispersal and transport of pollutants in two ways. 
Contaminants are first distributed fairly uniformly through- 
out the nearsurface air column, or marine layer. Further 
vertical transfer is then inhibited by the inversion, while 
lateral mixing continues until homogeneity is attained. 
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In winter, the above pattern is completely disrupted by 
the passage of migratory cyclones. However, during these 
times the marine layer is extremely deep and pollution is 
not nearly as much of a problem. 
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IV. MODEL DEVELOPMENT 



A. BASIC ASSUMPTIONS 

Since pollutants occur at concentrations of no more than 
a few parts per million in the Monterey Bay atmosphere, 
their presence can be assumed to have no effect on the mete- 
orological conditions. Additionally it must be recognized 
that the incorporation of spatial photochemistry (chemical 
reactivity of pollutant species) in a dispersion model is a 
complete subject for research in itself. Thus the concen- 
tration of a relatively nonreactive constituent, sulfur 
dioxide, is used in this study to model the diffusion of all 
contaminants. If the basic model can accurately predict SO^ 
concentrations, chemical constraints may later be considered. 

Atmospheric flow is assumed to be incompressible and to 
obey the laws of conservation of mass, momentum and energy. 
Wind is given as a function of x, y and t, but is considered 
to be independent of height and invariant in the vertical 
coordinate only above a friction layer, the depth of which 
is related to surface roughness. 

B. GAUSSIAN MODEL 

The concentration of SO 2 in the atmosphere must satisfy 

the continuity equation, written in flux form as: 

^ + V* (VC) = R + Q + DV^C (1.1) 

9t 
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where 



C = concentration of SO 2 at a point; 

R = time rate of change of concentration due to photo- 
chemical reaction; 

Q = source intensity; 

D = molecular diffusivity of SO 2 . 

The evolution of the velocity components u, v and w can be 
approximated by the Navier-Stokes equations for a nonrota- 
ting, incompressible Newtonian fluid. Conversion of these 
equations from Cartesian tensor form as expressed by Monin 
and Yaglom [1965, Eq. 1.6] to vector flux form yields: 

+ (V‘V)V = -gk - -Vp + vV^V . (1.2) 

8t P 

The system of equations is completed by the continuity 
equation : 

V*V = 0, (1.3) 

and the equation of state : 

p = pRT . (1.4) 

The wind velocity, V in equation (1.1), may be written 
as the sum of mean and eddy flow components: 

^ = V + V' (1.5) 

As used here, V is defined as the average wind computed over 
a volume characteristic of the mesh length to be used in the 
model. Since we assume sulfur dioxide to be photochemical ly 
inert, equation (1.1) may be written: 



The SO 2 concentration can also be considered the sum of 
deterministic (C) and stochastic (C*) components. Following 
the procedure delineated by Haltiner and Martin [1957], 
equation (1.6) is averaged to derive the equation for the 
mean concentration: 

+ V*VC + V*V'C = DV^C + Q (1.7) 

Note that while the average of a fluctuation alone or of the 
product of a fluctuation with a mean vanishes, the covariance 
of two fluctuations is not negligible. The covariance term 
is the turbulent flux, which is assumed proportional to the 
gradient in C (using the Prandtl Mixing Length principle; 
Haltiner and Martin, 1957, p. 245); 

VV'C = -V*KVC (1.8) 

where K is the turbulent diffusion coefficient. 

Observations have shown that diffusion due to molecular 
agitation is several orders of magnitude smaller than dif- 
fusion by turbulent eddies. Molecular diffusion can there- 
fore be neglected, and if the initial flow surface is 
relatively smooth, equation (1.7) reduces to: 

— + V* (VC) = 7*KVC + Q (1.9) 

at 

Several further approximations are required in formula- 
tion of the Gaussian model. The zonal mean velocity, u, is 
constant, and v=w=0. Since downstream advection is much 
greater than diffusion, diffusion along the plume axis is 
neglected (K =0) . K and K are assumed constant. For an 
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elevated point source of height H over a perfectly reflect- 
ing terrain (9C/9z = 0 at z=0 ) , the solution to equation 
(1.9), with the above assumptions, is; 



C= ^ exp(-y^/2aM {exp(-(z-H] ^/2a^ ) +exp (- [z+H] ^/2a ^) } 

— y z 2 



The quantities and are the standard deviations of 
particle distributions in the indicated directions and are 
functions of downwind distance and the statistical turbu- 
lence spectrum. 

Equation (1.10) expresses the concentration of SOj at 
ground level as a function of the height of the source. Due 
to the buoyancy and exit velocity of the stack gases, how- 
ever, the effective source height is really 



where is the stack height and AH is the rise of the plume 
after leaving the stack. 

A number of formulae have been developed to describe 
plume rise, most based upon empirical relationships. From 
observation of more than 700 stacks Moses and Carson [1968] 
developed the following expression: 




( 1 . 10 ) 



H = Hg + AH 



( 1 . 11 ) 




( 1 . 12 ) 



u 



u 



where A and B are constants which depend upon atmospheric 
stability. Plume rise figures used in this thesis are 
given in Table I. 



TABLE I 



PLUME RISE 

WIND VELOCITY ATMOSPHERIC STABILITY CLASS 





1 & 2 


3 


4 


5 


(M/SEC) 


(unstable) 


(slightly 

unstable) 


(neutral) 


(slightly 

stable) 


1. 


3206.51 


2237.42 


1479.85 


1090.52 


2. 


1603.25 


1118.71 


739.92 


545.26 


3. 


1068.86 


745.81 


493.28 


363.51 


4. 


801.63 


559.35 


369.96 


272.63 


5. 


641.30 


447.48 


295.97 


218.10 


6. 


534.42 


372.90 


246.64 


181.75 


7; 


458.07 


319.63 


211.41 


155.79 


8. 


400.81 


279.68 


184.98 


136.32 


9. 


356.28 


248.60 


164.43 


121.17 


10. 


320.65 


223.74 


147.98 


109.05 


11. 


291.50 


203.40 


134.53 


99.14 


12. 


267.21 


186.45 


123.32 


90.88 


13. 


246.65 


172.11 


113.83 


83.89 


14. 


229.04 


159.82 


105.70 


77.89 


15. 


213.77 


149.16 


98.66 


72.70 


16. 


200.41 


139.84 


92.49 


68.16 


17. 


188.62 


131.61 


87.05 


64.15 


18. 


178.14 


124.30 


82.21 


60.58 


19. 


168.76 


117.76 


77.89 


57.40 


20. 


160.33 


111.87 


73.99 


54.53 
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C. AN ALTERNATE MODEL 



As discussed earlier, the complex and irregular nature 
of turbulent boundary layer flow within the coastal airshed 
precludes the use of exact or measured meteorological data 
to describe diffusion. Nonetheless there exists the need 
to statistically parameterize this turbulence and include 
such considerations in pollutant dispersal models. The 
assvunption of homogeneous, isotropic flow which yields the 
Gaussian solution also eliminates flow boundaries from con- 
sideration. Wind velocity is required to be constant, a 
dangerous oversimplification. 

Obviously, the need exists for a multipurpose diffusion 
model capable of handling variable meteorological parame- 
ters. Since present knowledge is not adequate to permit a 
rigorous interpretation of turbulent transfer processes, the 
optimum approach should consider each aspect of the problem 
separately. The best approximation derived for each phase 
can then be incorporated in the model. 

Returning to equation (1.9), let horizontal diffusion be 

isotropic (K„=K, =K„) . The further assumption that K„ is 
X y n n 

constant in space and time is poor, since in fact the eddy 
diffusivity depends on both transport time and height; how- 
ever, this simplification can later be rectified. Thus 
equation (1.9) can be written 

“ = + Q (1.13) 

9t 9z 
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where the vertical advection, w9C , has been omitted as a 

3z 

first approximation (since w=0 near the ground except for 
terrain effects) . 

1. Boundary Conditions 

The boundary conditions are assumed to be: 

VC=0 at x=y=z=0, z=2 (1.14) 

where Z = inversion height. Exact lateral boundary condi- 
tions cannot be formulated, but as Shir and Shieh [1973] 
noted, the lack of well-posed boundary conditions does not 
cause serious problems. 

2 . Formulation of Numerical Model 

The accuracy of an analytical solution to equation 
(1.13) is restricted by the nonlinear nature of pollutant 
advection. In like manner, finite difference approximations 
to the advective term can introduce phase and amplitude 
errors in the numerical model. Simple upstream difference 
schemes, for example, produce a false dispersive effect 
which would make any results meaningless. Solutions obtain- 
ed using centered-difference forms of the governing equations 
are subject to stability limitations on the size of the time- 
step relative to the spatial mesh size of the model. These 
stability restrictions can lower computational efficiency by 
imposing a smaller timestep than would otherwise be desir- 
able. Thus a key disadvantage of explicit finite difference 
approximations to equation (1.13) is that the maximum time- 
step is fixed by the spatial grid distance rather than by 
the rate of change of the physical variables involved. 
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3 . Implicit Methods 

The unconditional stability of implicit methods 
makes them particularly desirable from a meteorological 
standpoint. Until recently, however, the implicit solution 
of large-matrix partial differential equations was accom- 
plished almost exclusively by iterative relaxation tech- ' 
niques. In pure form, such as with Liebmann Overrelaxation, 
these methods require large amounts of core storage and do 
not take advantage of the zero structure of solution mat- 
rices. 



methods known as alternating direction implicit (ADI) has 
been developed. Douglas [1962] modified the Crank-Nicholson 
equation to yield a three-step ADI scheme for the solution 
of mildly-nonlinear, three-dimensional parabolic and ellip- 
tic problems. The original alternating direction concept 
has since been expanded and generalized; a discussion of 
various techniques is given by Mitchell [1969] . 



To remove these objections, a noniterative class of 



Following Douglas' technique, the difference equa- 
tions for equation (1.13) are; 



X-SWEEP ; 




y-SWEEP : 



1 

At^2^y 




(1.15) 
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2-SWEEP : 



[ [iv^K V,]C^ 

At 2 2 ‘2 z z z" 



o 

As used here, the operators V and A apply at the point i, 
j, k. Space subscripting is eliminated to simplify notation. 
Forms of the finite-difference operators used are: 



u pH 

''x(y,z)^ijk 






n 



i+i / j / k~*^i-i , j ,k ^ ^ 
2Ax 9x 



n 



+ 0(Ax) 



i/ j/k 



(1.16) 



-.n 



^x(y,z)'“ijk 



,-2C?^i,+C? 



_ ^i+ 1 , jk~^’“i jk~*~^i- 1 , j , k _ 



(Ax) 



3^C 

3x^ 



n 



+ 0(Ax) 



if j#k 



When applied at successive grid points, equation (1.15) gen- 
erates a tridiagonal system of equations for The 

formal temporal accuracy of the solution is increased by the 
use of Crank-Nicholson centering about (i.e. by aver- 

aging terms on the RHS of equation (1.13) between the nth 
and (n+l)st levels). 
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V. MODEL APPLICATION 



A. SOURCE DATA 

The intensity of stack SO^ discharges and the physical 
design of the hypothetical power plant used in both thesis 
models are taken directly from data provided by Mr. R. C. 
Puckett of Pacific Gas and Electric Co. , Moss Landing. 

These data are summarized in Table II. 

B. GAUSSIAN MODEL 

1. Grid 

The southern half of Monterey Bay is divided into a 
receptor grid comprising 875 points: 25 points spaced one 
kilometer apart oriented east-west, by 35 points one km 
apart oriented north-south. The location of the power plant 
source is fixed at X=0.0, Y=0.0 and is the northwesternmost 
grid point. 

2 . Meteorological Data 
a. Wind Field 

For initial runs the surface wind direction is 
assumed constant at 315®. Wind velocity is varied from one 
mps to twenty mps in increments of one mps. 

Standard procedure, in estimating the diffusion 
of an effluent at a downwind distance of more than a few 
kilometers, is to assume a constant wind field. The center- 
line of the plume is then assumed to move in a straight line 
following the same trajectory computed at the source location. 
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TABLE II 



THEORETICAL POWER PLANT SPECIFICATIONS 

Location of Power Plant 
X-COORD y-COORD 

0.0 M 0.0 M 

Fuel: Oil (Primary) and Nat. Gas (as available) 

Actual Stack Height; 155.0 M 
Inside Stack Diameter; 6.0 M 
Stack Gas Exit Velocity; 20 M/SEC 
Stack Gas Temperature; 425 Deg K 
Stack Gas Constituents (in percent) ; 

NOj^; 72.1 
SO 2 ; 0.3 
COj; 13.7 
Particulates; xx.x 
Hydrocarbons ; xx.x 

Total Volume of Stack Gas Discharge; 7.01 x 10^ CU FT/HR 
Waste Heat Discharge; 1.0 x 10® BTU/HR 
SO 2 Source Intensity (est) ; 30 GM/SEC 
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This simple procedure is particularly inadequate in the 
coastal zone because of the complicated effects of the land- 
sea boundary and topographic influences. 

b. Atmospheric Stability and Dispersion Coefficients 

The standard deviations, and of plvime con- 

centration distribution in the horizontal and vertical vary 
with the turbulent structure of the marine layer. The dis- 
persion coefficient curves from Turner [1969] are used in 
the Gaussian model. Turner based his values of these param- 
eters upon stability estimates of the lowest few hiindred 
meters of the atmosphere, over "relatively open country". 

Stability is in turn determined by wind speed and insola- 
tion. Five classes are delineated (Table III) , with lapse / 

rates approaching neutral more closely for greater wind 
speeds . 

As previously shown, the thermal structure of 
the marine layer produces good vertical diffusion below the 
elevated inversion with virtually zero diffusion through the 
base of the stable layer. In a study of the same Moss Land- 
ing plant, Kraft [1971] noted the relative insensitivity of 
the atmosphere to the large quantities of waste heat dis- 
charged from the stacks. Thus the thesis model treats the 
inversion layer as a physical barrier to pollutant diffusion 
and ignores any local "heat island" effect of the plant. 
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TABLE III 

STABILITY CLASSIFICATION 



WIND SPEED 
(mps) 

<2 

2-3 

4-5 

6 

>6 



INSOLATION 



3 (strong) 



2 (moderate) 1 (weak) 



0 (night) 



12 2 
2 2 3 

2 3 3 

3 4 4 

3 4 4 

[after Turner, 1969] 



5 

5 

4 

4 

4 
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C. ALTERNATING DIRECTION IMPLICIT MODEL 
1. Grid and Timestep 

In order to provide for a more detailed analysis and 
to prevent eddy-scale phenomena from being smoothed, the 
horizontal mesh length for the ADI model was reduced to 0.5 
km. A 60-meter vertical spacing is used, dividing the ma- 
rine layer (600 m thick) into ten equal sublayers. The 
horizontal grid is oriented north-south (30 by 30 points) . 
The power plant is located at X=0.0, Y=7.5 km, Z=300 m. 

Selection of a timestep represents a critical deci- 
sion. The most important factor to be considered in deter- 
mining the step size is the time scale of atmospheric 
motions being described. Holland [1967] categorized the 
following time and length scales of turbulent atmospheric 
motion ; 

1) inertial turbulence: <10 sec, <50 m 

2) mechanical turbulence: 10-100 sec, 50-500 m 

3) thermal turbulence: 100-1000 sec, 500-5000 m. 

The mixing time for a typical turbulent marine layer 

has been described by Shieh (personal communication) as 
about ten minutes and by Lamb (personal communication) as 
15 minutes. Thus to parameterize the process of diffusion 
on a turbulent scale, a time increment of 30 seconds has 
been used. 
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2. Meteorological Conditions 
a. Wind Field 

For the ADI model, the wind vector 
Vij = u^j(x,y) + v^j(x,y) 

is required at all grid points for each sweep. As with the 
Gaussian method, a constant wind field is input, with all 
levels normalized to the input field (selected as 10 m) . 

For ease of computation the u-component of the wind is de- 
fined such that flow is along the axis of the "average" 
plume. Since total displacement over total travel time is 
a measure of the mean velocity, the displacement timescale 
is a natural one over which to define the average concentra- 
tion expressed by equation (1.13) . Thus the velocity com- 
ponents causing diffusion are deviations definable by 
statistical moments such as diffusion coefficients. 

Since wind data for specific levels within the 
marine layer are not available, the direction of horizontal 
wind flow is assumed constant with height. This assumption 
neglects the effects of directional shear; unfortunately the 
detail of available meteorological observations is not suf- 
ficient to correct this omission. 

Vertical wind profiles are assumed to be of power 

law form: 

where Vj^=V(10m) and Zj^=10m. The power coefficient, p, is 
varied from 0.1 to 0.4 according to stability class. Figure 
3 shows vertical wind profiles (normalized to the ten meter 
level) for each stability category. 
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NORMALIZED WIND VELOCEY U(Z)/U(10) 



FIGURE 3 ADI VELOCITY PROFILES 
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b. Atmospheric Stability 

In previous discussion the rate of diffusion of 
airborne pollutants has been attributed to the scale and 
intensity of atmospheric turbulence. This turbulence is 
directly related to stability. 

Classifications of atmospheric stability such as 
the Pasquill-Gif ford scheme used in the Gaussian model do 
not satisfactorily describe the topographic influences so 
important in generating turbulence. The Richardson nviraber, 
Ri and the Monin-Obukhov length, L are frequently-used sta- 
bility functions in boundary- layer theory. Unfortunately, 
both parameters depend upon the accuracy of input vertical 
temperature gradients. Recently, however, Golder [1972] 
related the Pasquill-Gif ford stability classes to the Monin- 
Obukhov length, incorporating surface roughness (Zq) and 
eliminating the requirement for vertical temperature pro- 
files. The Monin-Obukhov stability length is more conven- 
ient to use in the formulation of diffusion coefficients 
since it can be considered independent of altitude in the 
boundary layer. The Golder formula is 

(±) [0.216586-ln(1.2+— ) ] ^-10^^®^ (2.2) 

ij Z 0 

where f(s) = , s = stability class. 

(l+1.3ls| “• 

The sign of s determines the sign of the RHS of (2.2); nega- 
tive values are unstable, positive values stable. Golder 's 
formula is applied in the ADI model using a value of 1.5 
meters for z g . 
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c. Diffusion Coefficients 



Calculation of eddy dif fusivities is at best a 
guessing game, since data are sparse and estimates of both 
and vary in recent literature by at least an order of 
magnitude. These parameters are known to depend upon alti- 
tude, wind velocity, atmospheric stability, surface rough- 
ness and travel time but the relationship in each case is 
poorly established. 

Horizontal shearing stress has already been 
assumed constant. Travel time thus is probably the most im- 
portant factor in determining the coefficient of horizontal 
eddy diffusion, Kjj. As a pollutant cloud spreads downwind, 
the frequency of the turbulent elements most effective in 
dispersing particles decreases, so that the magnitude of 
eddy exchange increases. In order to parameterize this 
change, Kjj is initially set at lO^m^/sec (a frequently- 
quoted value) , and is then incremented within the time loop 
by lOOm^/sec at successive downstream grid points. When the 
value Kjj=10 ^m^/sec is reached incrementation stops; Kjj is 
held constant beyond that time. 

The vertical distribution of K under neutral 

z 

conditions has been given by Shir and Shieh [1973] as 

K = u*l (2.3) 

where 1 = koZ exp(-4 z/H) , 

u* = friction velocity, 
ko = von Karman constant = 0.41, 

H = depth of boundary layer (taken as 100m) . 
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(2.4) 



For non-neutral conditions, 

where the subscript s denotes values derived from the input 
"surface" wind velocity field, and the Monin-Obukhov length 
computed above. Luers [1973] describes the mean wind for 
the surface boundary layer (corresponding to 

u = ^[ln(^^) + (2.5) 

ko z L . 

from which u* can be computed. For neutral stability, s=0, 
i|;(z/L)=0, but for unstable conditions, s<0, 

z/L ^ , 

ip(^) = f d(f) (2.6) 

^ Zq/L^ ^ ^ 

For stable conditions, s>0, 

= 5.2(|-). 

Li JLi 

Using the value computed for u^, assuming the mean wind level 

is 10m, K -K (z=10m) is then calculated and extrapolated for 

(O-lOOm) . Above 100m the change in is assumed to be 

negligible and K =K (z=100m) . The normalized profiles for 

z z 

under neutral and slightly unstable conditions are dis- 

played in Figure 4. 
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FIGURE 4 ADI DIFFUSION COEFFICIENTS 
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VI. RESULTS 



A. GAUSSIAN MODEL 

Pollutant concentrations calculated by the Gaussian dif- 
fusion program under various meteorological conditions are 
described by Figures 5-9. Numerical values given are peak 
ground-level concentrations in micrograms per cu. meter. 

Figure 5 compares individual profiles for stability con- 
ditions ranging from neutral (s=4) to moderately unstable 
(s=2) . Wind speed and marine layer thickness for the three 
sets of curves are the same. Based upon an initial source 
intensity of 30 gm/sec, computed surface concentrations are 
within a fairly narrow range; from 2 . 3ygm-cu.m~^ for s=4 to 
11. 9ygm-cu.m“ ^ for s=2. However, downwind distance to the 
peak concentration varies by an order of magnitude, describ- 
ing the response of dispersion rates to increased thermal 
turbulence in the atmosphere. 

Figure 6 relates typical profiles of cases for which 
wind and stability are held constant, varying only the in- 
version height. Since the peak surface concentration as 
described by equation (1.10) is not related to layer thick- 
ness, the SO 2 ground-level maxima are identical, as is dis- 
tance to the peak value. Downstream from the point of 
maximum concentration the reduction of contaminant is much 
slower for the low-inversion case due to trapping in the 
marine layer. 
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FIGURE 5 DOWNWIND SO 2 CONCENTRATION 

AS A FUNCTION OF STABILITY 
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FIGURE 6 DOWNWIND S02 CONCENTRATION 
AS A FUNCTION OF MARINE LAYER DEPTH 



Variation of wind strength affects ground-level values 
as shown in Figure 7. Distance to peak concentration varies 
inversely with wind speed; however, the maximum intensity is 
greater for the intermediate wind speed (4 mps) than for 
either the 2mps or 16 mps cases. 

Figures 8 and 9 are composites of all cases studied. In 
Figure 8 , the maximum ground-level concentration is plotted 
as a function of wind speed and stability class. For winds 
stronger than 5 mps concentrations depend primarily upon 
stability; below 5 mps there is a pronounced increase in de- 
pendence upon wind speed, and a reversal in the trend of 
reduced peak concentrations with increasing distance from 
the source. For the stable cases (s=5) , peak surface values 
continued to increase to the edge of the grid but were sig- 
nificantly lower than values obtained for neutral or unstable 
conditions . 

Figure 9 shows downwind distance along the plume axis to 
the peak ground-level concentration, again as a function of 
stability and wind speed. Notice that the dependence upon 
stability class is almost total for stronger winds but that 
wind speed again becomes significant below about 5 mps. 

B. ADI MODEL 

Figures 10-13 show results obtained from the ADI model. 
Note that the basic purpose of these computations is to ex- 
amine the capabilities of the alternating direction method 
for the implicit solution to three-dimensional flow equations. 
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FIGURE 7 DOWNWIND S02 CONCENTRATION 
AS A FUNCTION OF WIND VELOCITY 
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*values at 25 km, concentration still increasing 



FIGURE 8 PEAK GROUND-LEVEL CONCENTRATION: ygm/m^ 
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STABILITY CLASS 



*values at 25 km, concentration still increasing 



FIGURE 9 DISTANCE TO PEAK GROUND-LEVEL CONCENTRATION 

FROM SOURCE; km. 
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Thus, concentration values are not to be compared analyti- 
cally with the Gaussian results to judge the accuracy of 
either method. 

In the 'course of developing the numerical scheme, solu- 
tions were developed for the case of Fickian diffusion, with 
K^=Ky=constant and K^=constant. A vertical, axial-plane 
profile describes these results in the marine layer (Figure 
10) . The solution is given for- a short-term (30 minute) 
calculation using a 30-second timestep. Concentrations fall 
off rapidly downstream and show strong vertical gradients. 
Defining the plume boundary as the surface at which the pol- 
lutant concentration is ten percent of the axial value on 
the same crosswind plane, the plume travels approximately 
five km downstream before ground contamination occurs. 

Figures 11-13 describe the effects of varying stability 
and wind velocity on dispersion, plotted on the axial verti- 
cal plane. These calculations include the horizontal and 
vertical diffusion coefficients discussed in MODEL INPUT. 
With reduced atmospheric stability, plume axial concentra- 
tion decreases more rapidly, but ground-level contamination 
occurs earlier. At a wind speed of 10 mps, the plume edge 
reaches the surface about 3.5 km downstream under slightly 
unstable conditions but travels nearly 11 km before ground- 
ing under neutral stability. Reducing wind speed slows the 
decrease in axial concentration but increases the distance 
to the surface contamination point. 
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C. DISCUSSION OF RESULTS 



The results of both models can be interpreted to yield 
some important generalizations concerning pollutant trans- 
port and diffusion in the atmosphere. The Gaussian model 
predicts statistical patterns of concentration, while the 
stochastic model deals with the probabilistic nature of dif- 
fusion. Both approaches have relevance in the atmosphere. 

From the patterns described by Figures 5-9 it is appar- 
ent that peak ground concentrations predicted for winds in 
excess of five mps show much less dependence upon the actual 
wind speed than do values obtained for winds less than five 
mps. This difference can be attributed in part to the pre- 
dominance of thermal turbulence at low wind speeds. Since 
the stability classification used in this study is derived 
as a function of wind velocity, results are difficult to 
interpret, but the fact seems well-established that weaker 
winds are more affected by thermal instability in the air. 
For strong winds the velocity profile must be dominated by 
the effects of mechanical turbulence, since only neutral or 
slightly unstable conditions are permitted at wind speeds in 
excess of five mps. This observation correlates well with 
the observations of Panofsky [1964] who found the high- 
frequency end of the turbulent velocity spectrum to be in- 
dependent of stability, and the low-frequency range strongly 
influenced by stability. In the case of low wind speeds, 
the simple power law formula used to define the wind field 
for the ADI model would therefore cease to be valid due to 
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the downwind dominance of thermally-driven turbulence in the 
lower atmosphere. 

In addition, under conditions of light wind the effective 
stack height is greatly increased by plume rise, frequently 
carrying the stack effluent through the inversion layer. 
Pollution hazards under such conditions are therefore not 
likely to be very serious. When an exceptionally strong in- 
version exists under stagnating anticyclonic flow, however, 
the potential for significant atmospheric contamination is 
much greater. 

Both methods predict more rapid ground-level contamina- 
tion as wind velocity increases. For lower wind speeds the 
Gaussian plume is released from a much greater effective 
height and the downwind dislocation of surface maxima is the 
expected result. However, the plume rise correction is not 
applied to the ADI model; thus the inverse relationship be- 
tween wind speed and ground contamination distance must be 
attributed to more effective vertical diffusion, accomplished 
by the increase in due to a larger u^ in equation (2.1) 
(see Figures 11 and 12) . 

If the size of the largest eddies responsible for 
spreading the plume vertically is considered to be approxi- 
mately the same as the crosswind plume dimensions one sees 
that mechanical turbulence can easily accomplish the neces- 
sary diffusion within the scale depth of the marine layer. 
Horizontal mixing is a much slower process, involving time- 
scales of hours or days until total equilibrium of concen- 
tration is established. Thus the turbulent response of the 
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boundary layer to topographic obstacles and other surface 
features will significantly alter vertical diffusion rates 
but will average out in horizontal dispersal. 

Concentration patterns produced by varying atmospheric 
stability (Figures 5, 8-9, 11-12) reveal the significance of 
vertical mixing determined by the marine layer thermal 
structure. Both models effected more rapid vertical dif- 
fusion in less stable air; the Gaussian model by means of 
larger values, the ADI model by means of increased turbu- 
lence as represented by the diffusion coefficient. These 
effects induce more rapid vertical diffusion, so that sur- 
face contamination at fairly high levels occurs within a few 
kilometers downstream of the source. Under stable conditions 
the damping of turbulence is so strong that flow is nearly 
laminar, trapping the pollutant in a duct near the effective 
height of injection and allowing essentially no pollutant 
ground deposition. Neutrally-stable marine layers yield 
results in between, with low- intensity surface pollution 
occurring at great distances downwind. 

Time restrictions precluded tuning the ADI model for 
study of marine layer depths other than 600 meters. How- 
ever, from the results of the Gaussian model (Figure 6) it 
is evident that significant pollution episodes are far less 
likely when the vertical extent of the marine layer is great. 
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VII. CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE STUDY 



The results of this thesis, though admittedly prelimi- 
nary, permit some broad conclusions regarding the general 
role of the atmosphere in dispersing contaminants. Non- 
reactive constituents and small particulates emitted into 
the atmosphere drift with the prevailing airstream, spread- 
ing progressively in both vertical and horizontal directions 
until they are removed by various natural processes (not 
considered in' this study) . The pollutants tend to be dis- 
tributed initially over a volume of air directly proportion- 
al to the wind speed, with simultaneous dilution by turbulent 
and convective motions which disturb the mean flow. Down- 
wind from a continuous emission, correlation of source data 
and various parameters of the idealized atmosphere shows 
peak concentrations to be directly proportional to rate of 
emission and inversely related to wind velocity and rate of 
diffusion. 

Effective emission height is shown to be a significant 
characteristic determining downstream concentrations. At 
greater elevations, pollutants must travel farther with 
correspondingly greater vertical and crosswind dilution be- 
fore ground-level contamination occurs. Since the effective 
source elevation is largely determined by wind velocity, the 
latter (along with inversion height, which is wind speed- 
related) is the most critical factor controlling pollution. 

In particular, under conditions of light wind when advection 
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is reduced and diffusion is uncertain, estimates of the dis- 
persive action of the atmosphere are not reliable. Strong 
winds produce concentration patterns which can be more de- 
pendably estimated from statistical fluid properties. 

The dependence of eddy diffusion rates upon atmospheric 
thermal stability is well-established and has been reaffirmed 
by results presented in this thesis. A first objective for 
further study should therefore be the definition of a sta- 
bility parameter which is a continuous function of marine 
layer stability. If accurate predictions are desired, eddy 
exchange coefficients must be computed from stability curves 
based on real-time, detailed thermal data (usually not rou- 
tinely available) . 

With specific regard to the downstream variation of eddy 
diffusion, literature is virtually nonexistent. The in- 
creased diffusive effectiveness of progressively larger tur- 
bulence elements acting upon a cloud or plume moving downwind 
seems obvious. The nature of such variance (treated linearly 
in this study) is not established and appears worthy of 
further evaluation. 

Any complete study of atmospheric diffusion in the 
coastal zone should consider the rapid changes in marine 
layer thermal structure occurring as a result of the temper- 
ature and terrain differences between land and water. Fu- 
ture research applied specifically to the Monterey area 
should study the effects of the land-sea breeze thermal cir- 
culation on pollutant dispersal. Under conditions of 
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onshore flow the thermal irregularities of underlying ter- 
rain will introduce serious complications into the analysis 
of diffusion; thus the effects of elevation, ground cover 
and moisture should be considered. 

Reasons for the preference of Gaussian dispersal models 
over more sophisticated schemes have already been discussed. 
The statistical method used in this study computes a solu- 
tion at 875 grid points for each of 20 wind speeds under 
four categories of insolation in less than eight seconds, 
requiring 49K bytes of core storage. Concentrations calcu- 
lated by the time-dependent ADI model require nearly 300 
seconds and more than 170K for each wind speed and stability 
class (9900 grid points) . Although no verifying data can be 
applied in either case, for short-term calculations the use 
of numerical simulation schemes is difficult to justify in 
light of the above cost comparisons. On the other hand, the 
terrific potential of time-dependent models to incorporate 
the effects of an observed, spatially- and temporally-vari- 
able wind field adjusted for topographic influences may out- 
weigh considerations of computational economy. 

The analysis of contaminant dispersal in the marine 
layer can provide valuable insight into the dynamics of that 
region. The history of an inert pollutant particle depends 
to a very great extent upon the turbulent properties of the 
near-surface atmospheric flow. Detailed wind observations 
and thermal profiles are rarely plentiful enough to reveal 
the strongly topography-dependent characteristics of the 
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turbulence theory. Observations of the transport of atmos- 
pheric contaminants such as sulfur dioxide may be used 
qualitatively to extend the available data base for the 
analysis of boundary-layer phenomena. 
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APPENDIX A 



ATMOSPHERIC SULFUR DIOXIDE; GENERAL CONSIDERATIONS 

A. SOURCES 

Coinbusion of sulfur-containing hydrocarbon fuels is the 
only significant cause of sulfur dioxide pollution in the 
atmosphere. Coal-burning power and industrial plants and 
the internal combustion engine are the major contributors. 
The governing reaction is 

S + O2 ^ SO2 . (A.l) 

B. ATMOSPHERIC REACTIVITY 

The oxidation of SO2 in a mixture of gases is slow, but 
the rate increases in sunlight and in the presence of NO2 , 
oxidants, and metallic oxides which act as catalysts. The 
reactions are: 



2SO2 t 02"^ 2SO3 


(A. 2 ) 


SO3 + H2O H2S0^ 


(A. 3 ) 


EFFECTS 




1 . Damage to Inorganic Materials 





Many substances suffer serious corrosive damage when 
exposed to sulfuric acid vapor, which is the secondary pol- 
lutant derived from SOj (equation A. 3 ). A few of these 
materials are listed below: 

Building stone: leaching, deterioration, discolora- 

tion; 
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Fabric, leather, dyes; weakening, embrittlement, 
fading; 

Metals; surface deterioration, net loss of metal 
(esp. zinc) ; 

Paper; fiber deterioration; 

Paint; discoloration; 

Rubber; cracking, loss of elasticity. 

2 . Damage to Vegetation 

A large number of plants, particularly some leafy 
commercial species, suffer serious leaf damage (bleaching 
and chlorosis) under the action of atmospheric SOj. The 
injury threshold varies according to species but for the 
more susceptible, such as pear trees, the maximum level for 
continuous eight-hour exposure is about 0.3 ppm. 

3. Effects on Animals 

As with plants, animals differ markedly in their 
susceptibility to inspired SOj. Primarily, SOj acts as an 
irritant to the respiratory system, damaging mucous linings 
of the upper tract. Relatively high levels (3.5 ppm) are 
necessary to produce a health hazard to the healthy adult 
human, but the aged (and those suffering from chronic re- 
spiratory or cardiovascular diseases) are much more suscep- 
tible. 

D. CONTROL OF SO 2 POLLUTION 

Suggested methods for the reduction and control of SO 2 
emissions include: 

use of low-sulfur fuels; 
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desulfurization of fuels before burning; 
stack gas scrubbing; 

use of higher stacks (at present the most effective 
and least expensive method available for reduc- 
ing ground-level concentrations) . 

All are partial answers at present, and the high cost and 
limited availability of such systems urges the development 
of alternate solutions. The availability of a workable dif- 
fusion model could permit the application of a regional 
"preferred location" scheme; power-generating loads would be 
shifted from plants in a region of high potential accumula- 
tion to those in areas of reduced potential, based upon 
prevailing meteorological conditions. 
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